Analysis for rape cases
library(dplyr) # A package for data manipulation
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(sf) # Simple feature for R
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(spdep) # Functions and tests for evaluating spatial patterns
## Loading required package: sp
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(tidyr) # Tools to create tidy data
library(INLA) # Integrated Nested Laplace Approximation package
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Loading required package: foreach
## Loading required package: parallel
## This is INLA_22.12.16 built 2022-12-23 13:43:44 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - To enable PARDISO sparse library; see inla.pardiso()
library(ggplot2) # A package for creating maps and graphs
library(viridis) # A package providing color palettes
## Loading required package: viridisLite
library(patchwork) # A package to compose plots
# For tables in RMarkdown
library(knitr)
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
# Packages used for visualization
library(RColorBrewer) # A package providing colour palettes for shading maps
# and other plots
library(tmap) # A package for static and interactive maps
library(ggplot2) # A package that implements the grammar of graphics, which is a term used to
# break up graphs into semantic components, such as geometries and layers.
library(mapview) # A package for interactive maps
library(cowplot) # Add-on to ggplot. It provides features that help us with creating
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:patchwork':
##
## align_plots
# publication-quality figures
load("~/Desktop/AA project/DS1_CrimeUttarPradesh/CrimeUttarPradesh.RData")
knitr::kable function:kable(data %>%
group_by(year) %>%
summarise(rape_observed = sum(rape), rape_expected=sum(e_rape)), booktabs = T, caption = "Rape cases in Uttar Pradesh by year") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
| year | rape_observed | rape_expected |
|---|---|---|
| 2001 | 1956 | 1524.186 |
| 2002 | 1413 | 1570.059 |
| 2003 | 908 | 1615.931 |
| 2004 | 1395 | 1661.804 |
| 2005 | 1214 | 1707.676 |
| 2006 | 1313 | 1753.549 |
| 2007 | 1646 | 1799.421 |
| 2008 | 1869 | 1845.293 |
| 2009 | 1758 | 1891.166 |
| 2010 | 1532 | 1937.038 |
| 2011 | 2040 | 1982.911 |
| 2012 | 1960 | 2028.783 |
| 2013 | 3047 | 2074.655 |
| 2014 | 3462 | 2120.528 |
st_read from sf package and call the object as
carto_up_sf.carto_up_sf <- st_as_sf(carto_up)
ggplot2 package.ggplot() +
geom_sf(data = carto_up_sf, color = "lightblue", fill = "white") +
coord_sf() + #axis limits and CRS
theme_bw() + # dark-on-light theme
theme(axis.title = element_text(size = 14),
axis.text = element_text(size = 12))
In order to use the same data structure for both the space-only model and later the space-time model, a new set of data is formed by aggregating both the observed and expected counts over time. A Poisson-log linear model is then fitted, assuming a BYM2 model for the random effects. Let each areal unit \(i\) be indexed by the integers \(1, 2,\dots,N\). The model is specified as:
\[ \begin{eqnarray} O_i &\sim & \text{Poisson}(\rho_i E_i)\\ \eta_i &= & \log \rho_i = b_0 + b_i\\ \boldsymbol{b} &= & \frac{1}{\sqrt{\tau_b}}(\sqrt{1-\phi}\boldsymbol{v}_{*} + \sqrt{\phi}\boldsymbol{u}_{*})\\ \end{eqnarray} \]
where \(\boldsymbol{v}_{*}\) and \(\boldsymbol{u}_{*}\) are standardised versions of \(\boldsymbol{u}\) and \(\boldsymbol{v}\). In detail, we specify the spatial random effect \(\boldsymbol{b}\) using a re-parameterisation of the Besag-York-Molli'e prior [@Besag1991], which is a convolution of an intrinsic conditional autoregressive (CAR) model and an independent and identically distributed Gaussian model, where
\(u_i\) is the spatially structured component defined by an intrinsic CAR prior: \(\boldsymbol{u}\sim ICAR(\boldsymbol{W}, \sigma^2_u)\),
\(v_i\) the unstructured component defined with prior: \(v_s \overset{iid}{\sim} \text{Normal}(0,\sigma^2_v)\).
carto_up_sf_nb = poly2nb(carto_up_sf, queen=TRUE)
summary(carto_up_sf_nb)
## Neighbour list object:
## Number of regions: 70
## Number of nonzero links: 324
## Percentage nonzero weights: 6.612245
## Average number of links: 4.628571
## 2 regions with no links:
## 4 11
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8 9
## 2 2 2 13 13 16 12 6 3 1
## 2 least connected regions:
## 45 60 with 1 link
## 1 most connected region:
## 16 with 9 links
inla format using the
function nb2WB(). call the graph as
carto_india_sf.adjnb2INLA("carto_up_sf.graph", carto_up_sf_nb)
carto_up_sf.adj=paste(getwd(), "/carto_up_sf.graph", sep ="")
ID_area.data_agg = data %>% group_by (ID_area) %>%
summarize(rape_observed = sum(rape),
rape_expected = sum(e_rape)) %>%
dplyr::rename(rape_obs = rape_observed, rape_exp = rape_expected)
data_agg = data_agg %>% mutate (rape_SMR = rape_obs/rape_exp)
Remember that, before to produce the map, you need to join the sf object and the data frame data_agg. To do so you can use the function left_join from the library dplyr.
data_agg$rape_SMRcat = cut(data_agg$rape_SMR,
breaks=c(min(data_agg$rape_SMR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(data_agg$rape_SMR)), include.lowest = T)
map_SMR = left_join(carto_up_sf, data_agg, by = c("ID_area" = "ID_area"))
and plot:
ggplot() + geom_sf(data = map_SMR, col = NA) + aes(fill = rape_SMRcat) +
theme_bw() + scale_fill_viridis_d() +
guides(fill=guide_legend(title="SMR for rape cases"))
Map of the average rape SMRs over the period 2001-2014 in Uttar Pradesh
INLA.
Remember to create an ID for the areas, i.e. index vector, and use a
bym2 prior for the spatial random effect setting the PC
priors [@Simpson2017] seen in Session 2.2.
Monitor also the Watanabe-Akaike information criterion (WAIC) including
control.compute=list(waic=TRUE)ID = seq(1,70)
formula_BYM2 = rape_obs ~ f(ID, model="bym2", graph="carto_up_sf.graph",
hyper=list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3))))
sBYM.model = inla(formula=formula_BYM2, family="poisson", data=data_agg, E=E, control.compute=list(waic=TRUE))
## Warning in set.warn("E", nm): Argument 'E=E' expanded to NULL or gave an error.
## This might be an error and you are requested to check this out.
## Move on with default values...
## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'bym2' in section 'latent' is marked as 'experimental'; changes may appear at any time.
## Use this model with extra care!!! Further warnings are disabled.
#Relative risks
RR_sBYM = c()
for(i in 1:70){
RR_sBYM[i] = inla.emarginal(function(x) exp(x),
sBYM.model$marginals.random$ID[[i]])
}
#Posterior probabilities
RR_sBYM_marg = sBYM.model$marginals.random$ID[1:70]
PP_sBYM = lapply(RR_sBYM_marg, function(x) {1-inla.pmarginal(0,x)})
resRR_PP = data.frame(resRR=RR_sBYM,
PP=unlist(PP_sBYM),
ID_area=data_agg[,1])
ggplot2 package, produce a map of the posterior
mean of the residual RRs and the posterior probabilities that the
residual RRs are > 1resRR_PP$resRRcat = cut(resRR_PP$resRR, breaks=c(min(resRR_PP$resRR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(resRR_PP$resRR)),include.lowest = T)
PP that the residual
RRs is > 1, use the following breakpoints [0-0.2], (0.2-0.8],
(0.8-1].# breakpoints
resRR_PP$PPcat = cut(resRR_PP$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
sf object and data frame with the
posterior estimatesmap_RR_PP = left_join(carto_up_sf, resRR_PP, by = c("ID_area" = "ID_area"))
PP using ggplot2
package.p1 = ggplot() + geom_sf(data = map_RR_PP) + aes(fill = resRRcat) +
theme_bw() + scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatial model for rape cases") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
p2 = ggplot() + geom_sf(data = map_RR_PP) + aes(fill = PPcat) +
theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatial model for rape cases") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
p1|p2
As the BYM2 has the structured (CAR) and unstructured (iid) components it might be useful to get some ideas about the strength of the spatially structured components as this would indicate the level of clustering in the data. To do so we ca simply obtain the posterior summary of the phi hyperparameter. How do you interpret it?
sBYM.model$summary.hyperpar
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID 2.9275149 0.6013841 1.9151889 2.8700128 4.2739890 2.7605712
## Phi for ID 0.6384346 0.1726621 0.2743392 0.6548827 0.9171358 0.7153379
ANSWER: This tells us that about 2/3 of the spatial variability is explained by the spatially structured component - which makes sense if we look at the map of the RR, which show a degree of spatial clustering.
Now, we extend the above analysis to a separable space-time model without interactions. For the temporal component, we use the specification introduces in the lecture with a temporal unstructured random effect and a structured one (RW1 prior).
Let each areal unit \(i\) be indexed
by the integers \(1, 2,\dots,N\). As in
the spatial case, we use a Poisson distribution to model the number of
hospital admission \(O_it\), in area
\(i\) at time \(t\). The mathematical specification of the
model includes now an additional temporal dependence term, which can be
modeled using a non-stationary random walk prior: ${i} ({i-1},
^2_{}).
The model implement in this practical assumes no space-time interaction
and a spatial convolution with random walk in time:
\[ \begin{eqnarray} O_{it} &\sim & \text{Poisson}(\rho_{it} E_{it}) \\ \log \rho_{it} &= & b_0 + b_i + \gamma_t + \psi_t \\ \boldsymbol{b} &= & \frac{1}{\sqrt{\tau_b}}(\sqrt{1-\phi}\boldsymbol{v}_{*} + \sqrt{\phi}\boldsymbol{u}_{*})\\ \gamma_t & \sim & \hbox{RW(1)}\\ \psi_t & \sim & N(0,\sigma^2_{\psi})\\ \end{eqnarray} \]
As seen before \(\boldsymbol{v}_{*}\) and \(\boldsymbol{u}_{*}\) are standardised versions of \(\boldsymbol{u}\) and \(\boldsymbol{v}\). In detail, we specify the spatial random effect \(\boldsymbol{b}\) using a re-parameterisation of the Besag-York-Molli'e prior, which is a convolution of an intrinsic conditional autoregressive (CAR) model and an independent and identically distributed Gaussian model, where
\(u_i\) is the spatially structured component defined by an intrinsic CAR prior: \(\boldsymbol{u}\sim ICAR(\boldsymbol{W}, \sigma^2_u)\),
\(v_i\) the unstructured component defined with prior: \(v_s \overset{iid}{\sim} \text{Normal}(0,\sigma^2_v)\).
#Join the data with the shapefile so the order of the shapefile is maintained.
data_ST = left_join(carto_up_sf, data, by="ID_area")
#Rename the columns of Observed and Expected as we did before
data_ST = data_ST %>% dplyr::rename(rape_obs = rape, rape_exp = e_rape)
#Create the ID for year (time)
data_ST$ID.time = data_ST$year - 2001
#Create the ID for space
data_ST$ID.space = rep(seq(1,70),each=14)
INLA, monitoring the WAIC and using PC
priors (as specified in Session 2.2) for the bym2 and
rw1 models (note: for the rw1 model you could
try the following specification
f(ID.time,model="rw1", hyper=list(prec = list(prior = "pc.prec",param = c(0.5 / 0.31, 0.01))))).
Call the output as stBYM.modelformula_ST_noint = rape_obs ~ f(ID.space, model="bym2", graph=carto_up_sf.adj,
hyper=list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3)))) + f(ID.time, model="rw1", hyper=list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01))))
stBYM.model = inla(formula=formula_ST_noint, family="poisson", data=data_ST, E=E,
control.compute=list(waic=TRUE))
## Warning in set.warn("E", nm): Argument 'E=E' expanded to NULL or gave an error.
## This might be an error and you are requested to check this out.
## Move on with default values...
#Spatial Relative risks
RR_stBYM = c()
for(i in 1:70){
RR_stBYM[i] = inla.emarginal(function(x) exp(x),
stBYM.model$marginals.random$ID.space[[i]])
}
#Posterior probabilities (for spatial RR)
RR_stBYM_marg = stBYM.model$marginals.random$ID.space[1:70]
PP_stBYM = lapply(RR_stBYM_marg, function(x) {1-inla.pmarginal(0,x)})
#Temporal Relative risks and CI95
RR_stRW_RR = c()
RR_stRW_lo = c()
RR_stRW_hi = c()
for(i in 1:14){
#Posterior mean
RR_stRW_RR[i] = inla.emarginal(function(x) exp(x),
stBYM.model$marginals.random$ID.time[[i]])
#2.5% quantile
RR_stRW_lo[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
#97.5% quantile
RR_stRW_hi[i] = inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
}
RR_stRW = data.frame(RR=RR_stRW_RR,low=RR_stRW_lo,high=RR_stRW_hi)
RR_stWR)Temp1 = ggplot(RR_stRW, aes(seq(2001,2014), RR)) + geom_line() + ggtitle("ST model No Int") + geom_ribbon(aes(ymin=low,ymax=high), alpha=0.2) + labs(x="year")
Temp1
resRR_PP_st = data.frame(resRR=RR_stBYM,
PP=unlist(PP_stBYM),
ID_area=data_agg[,1])
# breakpoints
resRR_PP_st$resRRcat = cut(resRR_PP_st$resRR, breaks=c(min(resRR_PP_st$resRR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(resRR_PP_st$resRR)),include.lowest = T)
resRR_PP_st$PPcat = cut(resRR_PP_st$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
map_RR_ST = left_join(carto_up_sf, resRR_PP_st, by = c("ID_area" = "ID_area"))
p3 = ggplot() + geom_sf(data = map_RR_ST) + aes(fill = resRRcat) +
theme_bw() + scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title="RR")) + ggtitle("RR ST model") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
p4 = ggplot() + geom_sf(data = map_RR_ST) + aes(fill = PPcat) +
theme_bw() +
scale_fill_viridis(
option = "plasma",
name = "PP ST model",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP ST model") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
(p1|p2) / (p3|p4)
Spatio-temporal model: Map of the residual RRs and posterior probabilities for rape cases
Now, we extend the above spatio-temporal analysis to a separable space-time model with type I interaction:
\[ \begin{eqnarray} O_{it} &\sim & \text{Poisson}(\rho_{it} E_{it}) \\ \log \rho_{it} &= & b_0 + b_i + \gamma_t + \psi_t + \delta_{it} \\ \boldsymbol{b} &= & \frac{1}{\sqrt{\tau_b}}(\sqrt{1-\phi}\boldsymbol{v}_{*} + \sqrt{\phi}\boldsymbol{u}_{*})\\ \gamma_t & \sim & \hbox{RW(1)}\\ \psi_t & \sim & N(0,\sigma^2_{\psi})\\ \delta_{it} & \sim & \hbox{Normal}(0, \sigma^2_{\delta}) \end{eqnarray} \]
INLA. Remember
that you need to create an index which goes from 1 to the length of the
dataset (i.e. the space x time). For the iid model defining
the interaction term, use the PC prior previously used for the
rw1 model. Call the output stIntI.BYM.model
(remember to monitor the WAIC)data_ST$ID.space.time = seq(1,dim(data_ST)[1])
formula_ST_intI = rape_obs ~ f(ID.space, model="bym2", graph=carto_up_sf.adj,
hyper=list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01)),
phi = list(
prior = "pc",
param = c(0.5, 2 / 3)))) +
f(ID.time,model="rw1", hyper=list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01))))+
f(ID.space.time,model="iid", hyper=list(prec = list(
prior = "pc.prec",
param = c(0.5 / 0.31, 0.01))))
stIntI.BYM.model = inla(formula=formula_ST_intI, family="poisson", data=data_ST, E=E,
control.compute=list(dic=TRUE, waic=TRUE))
## Warning in set.warn("E", nm): Argument 'E=E' expanded to NULL or gave an error.
## This might be an error and you are requested to check this out.
## Move on with default values...
#Spatial Relative risks
RR_stIntI.BYM = c()
for(i in 1:70){
RR_stIntI.BYM[i] = inla.emarginal(function(x) exp(x),
stIntI.BYM.model$marginals.random$ID.space[[i]])
}
#Posterior probabilities (for spatial RR)
RR_stIntI.BYM_marg = stIntI.BYM.model$marginals.random$ID.space[1:70]
PP_stIntI.BYM = lapply(RR_stIntI.BYM_marg, function(x) {1-inla.pmarginal(0,x)})
#Temporal Relative risks and CI95
RR_stIntI.RW_RR = c()
RR_stIntI.RW_lo = c()
RR_stIntI.RW_hi = c()
for(i in 1:14){
#Posterior mean
RR_stIntI.RW_RR[i] = inla.emarginal(function(x) exp(x),
stIntI.BYM.model$marginals.random$ID.time[[i]])
#2.5% quantile
RR_stIntI.RW_lo[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), stIntI.BYM.model$marginals.random$ID.time[[i]]))
#97.5% quantile
RR_stIntI.RW_hi[i] = inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stIntI.BYM.model$marginals.random$ID.time[[i]]))
}
RR_stIntI.RW = data.frame(RR=RR_stIntI.RW_RR,low=RR_stIntI.RW_lo,high=RR_stIntI.RW_hi)
RR_stWR)Temp2 = ggplot(RR_stIntI.RW, aes(seq(2001,2014), RR)) + geom_line() + ggtitle("ST model Int I") + geom_ribbon(aes(ymin=low,ymax=high), alpha=0.2) + labs(x="year")
Temp1 | Temp2
RR_stIntI.BYM) with
ggplot2 package using the following breakpoints [min,0.4],
(0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max].
Compare this map against the map of the residual RR obtained from the
spatio-temporal model with no interaction. Try to comment the
resultsresRR_PP_stIntI = data.frame(resRR=RR_stIntI.BYM,
PP=unlist(PP_stIntI.BYM),
ID_area=data_agg[,1])
# breakpoints
resRR_PP_stIntI$resRRcat = cut(resRR_PP_stIntI$resRR, breaks=c(min(resRR_PP_stIntI$resRR),
0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
max(resRR_PP_stIntI$resRR)),include.lowest = T)
resRR_PP_stIntI$PPcat = cut(resRR_PP_stIntI$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
map_RR_ST.IntI = left_join(carto_up_sf, resRR_PP_stIntI, by = c("ID_area" = "ID_area"))
p5 = ggplot() + geom_sf(data = map_RR_ST.IntI) + aes(fill = resRRcat) +
theme_bw() + scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title="RR")) + ggtitle("RR ST model Int I") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
p6 = ggplot() + geom_sf(data = map_RR_ST.IntI) + aes(fill = PPcat) +
theme_bw() +
scale_fill_viridis(
option = "plasma",
name = "PP ST model Int I",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP ST model Int I") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
(p1|p2) / (p3|p4) / (p5|p6)
Spatio-temporal model: Map of the residual RRs and posterior probabilities for rape cases
ANSWER: We basically see that across the different models there is tiny difference in the spatial residuals. Let’s now look at the space-time interaction.
data_ST$intI = stIntI.BYM.model$summary.random$ID.space.time$mean
data_ST$intI_cat = cut(data_ST$intI, breaks=c(-1,-0.05,
-0.01, 0.01, 0.05, 1),include.lowest = T)
ggplot() +
geom_sf(data = data_ST, aes(fill = intI_cat))+ theme_bw() + scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title=NULL)) +
theme(text = element_text(size=20),
axis.text.x = element_blank(),
axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 3, labeller=labeller(ID.year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005", "6"="2006", "7"="2007", "8"="2008", "9"="2009", "10"="2010", "11"="2011", "12"="2012", "13"="2013", "14"="2014"))) +
labs("")
We can see that there is not clear pattern in the interactions.
dat.hyper2 =
round(
data.frame(median = stIntI.BYM.model$summary.hyperpar[,4],
LL = stIntI.BYM.model$summary.hyperpar[,3],
UL = stIntI.BYM.model$summary.hyperpar[,5]),
digits = 3)
row.names(dat.hyper2) =
rownames(stIntI.BYM.model$summary.hyperpar)
knitr::kable(dat.hyper2, caption = "Posterior median and 95% CrI of hyperparameters for rape cases.") %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
| median | LL | UL | |
|---|---|---|---|
| Precision for ID.space | 2.763 | 1.822 | 4.195 |
| Phi for ID.space | 0.723 | 0.302 | 0.942 |
| Precision for ID.time | 13.107 | 5.660 | 29.793 |
| Precision for ID.space.time | 12.289 | 10.481 | 14.173 |
dat.WAIC = data.frame(model = c("Spatial", "SpatTemp no int", "SpatTemp typeI"),
WAIC = round(c(sBYM.model$waic$waic, stBYM.model$waic$waic, stIntI.BYM.model$waic$waic))
)
row.names(dat.WAIC) = NULL
knitr::kable(dat.WAIC, caption = "WAIC of the fifferent models for rape cases") %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
| model | WAIC |
|---|---|
| Spatial | 649 |
| SpatTemp no int | 8092 |
| SpatTemp typeI | 6245 |
ANSWER: Although the spatial model has a lower WAIC it cannot be directly compared to the two spatio-temporal models, as it is essentially based on a different data set (i.e. data were aggregated removing the temporal component). The fair comparison in between the two spatio-temporal models, and we observe that the model with type I interaction performs better that the model without interaction.